

yearstring = [ '1998';'1999';'2000';'2001';'2002';'2003';'2004';'2005'];

st = 1;

en = 1;

for y = 1:length(yearstring),
    
    
    
    S = dir(['D:\Manoj\definitive_Don_41GB!\' yearstring(y,:) '\*.min']);
    
    for i= 1:length(S),
        
        filestring(i,:) = S(i).name;
        
    end;
    
    [b,i,j] = unique(filestring(:,1:3),'rows');
    
    en = st + length(b)-1;
    
    uniquestations(st:en,:) = b;
    
    st = en + 1;
    
    fprintf('Year = %s stations = %d data days = %d\n', yearstring(y,:),length(b),length(S)-2);
    
end;



%%



% [b1,i1,j1] = unique(uniquestations,'rows');

% there are 106 unique stations

%save c:\manoj\projects\gary\1_min_station_code b1;

load c:\manoj\projects\gary\1_min_station_code b1;

date_array = zeros([sum(daysinyear([1998:2005])),length(b1)]);

daysaxis = datenum(1998,1,1):datenum(2005,12,31);





for y = 1:length(yearstring),
    
    
    
    S = dir(['D:\Manoj\definitive_Don_41GB!\' yearstring(y,:) '\*.min']);
    
    for i= 3:length(S),
        
        filestring(i-2,:) = S(i).name;
        
    end;
    
    
    
    [b, i, j] = unique(filestring(:,1:3),'rows');%Find the unique stations
    
    [d, f, g] = intersect(b1,b,'rows');%Find the intersection with the global stations database
    
    thisdayaxis = datenum(filestring(:,4:11),'yyyymmdd');%Find the time axis for all the stations in this directory
    
    
    
    for ii = 1:length(f), %For length of number of unique stations
        
        L = j == ii; %Select the current station
        
        [c,ia,ib] = intersect(daysaxis, thisdayaxis(L));%Time axis for the current station
        
        date_array(ia,f(ii)) = 1; %Put all the available data time to one
        
    end;
    
    
    
    fprintf('Year = %s stations = %d data days = %d\n', yearstring(y,:),length(b),length(S)-2);
    
end;







%date_array gives a matrix of available stations





%%

more off;



load c:\manoj\projects\gary\1_min_station_code b1;

yearstring = [ '1998';'1999';'2000';'2001';'2002';'2003';'2004';'2005'];

a = sum(daysinyear([1995:2006]))*24*60;

start_date = datenum(1995,1,1);



for i = 3:length(b1),
    
    
    
    nc_fname = ['C:\Manoj\projects\gary\minute_data\' b1(i,:) '.nc'];
    
    S1 = dir(nc_fname);
    
    
    
    if isempty(S1),
        
        ncid = netcdf.create(nc_fname,'NC_NOCLOBBER');
        
        data_time_axis_dim = netcdf.defDim(ncid,'data_time_axis_dim',a);
        
        X_ID = netcdf.defVar(ncid,'X','int', data_time_axis_dim);
        
        Y_ID = netcdf.defVar(ncid,'Y','int', data_time_axis_dim);
        
        Z_ID = netcdf.defVar(ncid,'Z','int', data_time_axis_dim);
        
        netcdf.putAtt(ncid,X_ID,'_FillValue',int32(999999));
        
        netcdf.putAtt(ncid,Y_ID,'_FillValue',int32(999999));
        
        netcdf.putAtt(ncid,Z_ID,'_FillValue',int32(999999));
        
        netcdf.endDef(ncid);
        
    else
        
        ncid = netcdf.open(nc_fname,'NC_WRITE');
        
        X_ID = netcdf.inqVarID(ncid,'X');
        
        Y_ID = netcdf.inqVarID(ncid,'Y');
        
        Z_ID = netcdf.inqVarID(ncid,'Z');
        
    end;
    
    
    
    
    
    
    
    for j = 1:length(yearstring),
        
        
        
        S = dir(['D:\Manoj\definitive_Don_41GB!\' yearstring(j,:) '\' b1(i,:) '*.min']);
        
        
        
        if ~isempty(S),
            
            
            
            for k = 1 : length(S),
                
                fname = ['D:\Manoj\definitive_Don_41GB!\' yearstring(j,:) '\' S(k).name];
                
                data_struct = intermagnet(fname, 'XYZF');
                
                if(strcmp(data_struct.CORD,'XYZF')),
                    
                    display('HDZF to XYZF conversion');
                    
                end;
                
                fprintf('file name: %s\n',S(k).name);
                
                where_to_store = ( datenum(S(k).name(:,4:11),'yyyymmdd') - start_date ) * 1440 ;
                
                %netcdf.putVar(ncid,varid,start,count,data)
                
                % multiplied by 10 to preserve the decimal points
                
                netcdf.putVar(ncid, X_ID, where_to_store, 1440, data_struct.DATA(:,1) .* 10 );
                
                netcdf.putVar(ncid, Y_ID, where_to_store, 1440, data_struct.DATA(:,2) .* 10);
                
                netcdf.putVar(ncid, Z_ID, where_to_store, 1440, data_struct.DATA(:,3) .* 10);
                
                
                
                
                
            end;
            
        end;
        
    end;
    
    
    
    netcdf.close(ncid);
    
    
    
end;



%%

% reading 1 min binaru imag data



yearstring = [ '1995';'1996';'1997';'1998';'1999';'2000';'2001';'2002';'2003';'2004';'2005';'2006'];



a = sum(daysinyear([1995:2006]))*24*60;

start_date = datenum(1995,1,1);

fid = fopen('c:\manoj\temp\Gary_1_min_log_error.txt','at');



%S = textread('D:\Manoj\imag19922004\file_listing.txt','%s');



%S = textread('D:\Manoj\imag19922004\file_listing0506.txt','%s');

%S = textread('D:\Manoj\imag19922004\IMAG2005\missig_files.txt','%s');

%S  = dir('D:\Manoj\imag19922004\IMAGfromftp\*.zip');

S = textread('D:\Manoj\imag19922004\IMAGfromDVD\allfiles_edited.txt','%s');



% for i = 1:length(S),

%     filestring(i,:) = S(i).name;

% end;



for i= 1:length(S),
    
    filestring(i,:) = upper(cell2str(S(i)));
    
end;



filestring = filestring(:,1:8);%there may be zip and bin files with for the sam estation



[filestring, ia, ib] = unique(filestring,'rows');%Find the unique stations




for i = 93:length(b),
    
    
    
    nc_fname = ['C:\Manoj\projects\gary\minute_data\' b(i,:) '.nc'];
    
    S1 = dir(nc_fname);
    
    
    
    if isempty(S1),
        
        ncid = netcdf.create(nc_fname,'NC_NOCLOBBER');
        
        data_time_axis_dim = netcdf.defDim(ncid,'data_time_axis_dim',a);
        
        X_ID = netcdf.defVar(ncid,'X','int', data_time_axis_dim);
        
        Y_ID = netcdf.defVar(ncid,'Y','int', data_time_axis_dim);
        
        Z_ID = netcdf.defVar(ncid,'Z','int', data_time_axis_dim);
        
        netcdf.putAtt(ncid,X_ID,'_FillValue',int32(999999));
        
        netcdf.putAtt(ncid,Y_ID,'_FillValue',int32(999999));
        
        netcdf.putAtt(ncid,Z_ID,'_FillValue',int32(999999));
        
        netcdf.endDef(ncid);
        
    else
        
        ncid = netcdf.open(nc_fname,'NC_WRITE');
        
        X_ID = netcdf.inqVarID(ncid,'X');
        
        Y_ID = netcdf.inqVarID(ncid,'Y');
        
        Z_ID = netcdf.inqVarID(ncid,'Z');
        
    end;
    
    
    
    L = ib == i; %select a station
    
    filesubset = filestring(L,:);
    
    
    
    
    
    for j = 1:sum(L),
        
        
        
        
        
        binfname = ['D:\Manoj\imag19922004\IMAGfromDVD\MAG' filesubset(j,4:5) '\' filesubset(j,1:3) '\' filesubset(j,1:8) '.bin'];
        
        %binfname = ['D:\Manoj\imag19922004\IMAG2005\2005\' filesubset(j,1:3) '\' filesubset(j,1:9) 'bin'];
        
        
        
        %binfname = ['D:\Manoj\imag19922004\IMAGfromftp\' filesubset(j,1:9) 'bin'];
        
        
        
        S = dir(binfname);
        
        
        
        if isempty(S),
            
            zipfilename = ['D:\Manoj\imag19922004\IMAGfromDVD\MAG' filesubset(j,4:5) '\' filesubset(j,1:3) '\' filesubset(j,1:8) '.zip'];
            
            %zipfilename = ['D:\Manoj\imag19922004\IMAGfromftp\' filesubset(j,:)];
            
            unzip(zipfilename, ['D:\Manoj\imag19922004\IMAGfromDVD\MAG' filesubset(j,4:5) '\' filesubset(j,1:3) '\']);
            
            %unzip(zipfilename, ['D:\Manoj\imag19922004\IMAGfromftp\']);
            
        end;
        
        
        
        
        
        
        
        data_struct = INTERMAG(binfname, 'XYZF');
        
        
        
        if ~isempty(data_struct) && (max(abs(diff(data_struct.MINUTELY(:,1)))) < 1.1/(24*60)),
            
            
            
            data_struct.MINUTELY(isnan(data_struct.MINUTELY)) = 99999.9;    %no data (this will get converted to 999999 by int32(99999.9*10)
            
            fprintf('file name: %s\n',binfname);
            
            where_to_store = ( data_struct.MINUTELY(1,1) - start_date ) * 1440 ;
            
            length_of_record = length(data_struct.MINUTELY(:,1));
            
            
            
            %netcdf.putVar(ncid,varid,start,count,data)
            
            %multiplied by 10 to preserve the decimal points
            
            netcdf.putVar(ncid, X_ID, where_to_store, length_of_record,  int32 (data_struct.MINUTELY(:,2) .* 10) );
            
            netcdf.putVar(ncid, Y_ID, where_to_store, length_of_record,  int32 (data_struct.MINUTELY(:,3) .* 10) );
            
            netcdf.putVar(ncid, Z_ID, where_to_store, length_of_record,  int32 (data_struct.MINUTELY(:,4) .* 10) );
            
            fprintf(fid,'%s             OK \n',binfname);
            
        else
            
            fprintf(fid,'%s             ERROR \n',binfname);
            
            
            
        end;
        
        
        
    end;
    
    
    
    netcdf.close(ncid);
    
    
    
end;

fclose(fid);



%%

% Compare the WDC Edn 1 min data and imag data I have

%%

%compare list wdc edin & my data base



fid = fopen('c:\manoj\projects\gary\list_1min_wdc_edi.txt','rt');

fid1 = fopen('c:\manoj\projects\gary\list_1min_wdc_edi_edited.txt','wt');

station_name = 'KKK';

while ~feof(fid),
    
    S = fgetl(fid);
    
    if any(ismember([1995:2006],str2num(S(1:4)))),
        
        fprintf(fid1,'%s %s\n',station_name, S);
        
    elseif ~any(ismember([1950:2010],str2num(S(1:4)))),
        
        a = find(S == '(');
        
        if ~isempty(a),
            
            station_name = S(a(end)+1:a(end)+3);
            
        else,
            
            display('error');
            
            break;
            
        end;
        
    else
        
        
        
    end;
    
end;

fclose all;



%%



fid = fopen('c:\manoj\projects\gary\list_1min_wdc_edi_edited.txt','rt')

i=1;

while ~feof(fid),
    
    SS(i,:) = fscanf(fid,'%s',4);
    
    i=i+1;
    
end;

fclose(fid);

[b_edin, ia_edi, ib_edin] = unique(upper(SS(:,1:3)),'rows');%Find the unique stations



i_count = 1;

mon = ['JFMAMJJASOND'];

justmonths = ['JAN';'FEB';'MAR';'APR';'MAY';'JUN';'JUL';'AUG';'SEP';'OCT';'NOV';'DEC']





for i = 1:length(SS),
    
    for j = 1:12,
        
        if SS(i,7+j) == mon(j),
            
            new_S(i_count,:) = [SS(i,1:3) SS(i,6:7) justmonths(j,:)];
            
            i_count = i_count + 1;
            
        end;
        
    end;
    
end;





%%

fid = fopen('D:\Manoj\imag19922004\IMAGfromDVD\allfile_edited.txt','rt')

i=1;



while ~feof(fid),
    
    S(i,:) = fscanf(fid,'%s',1)';
    
    i=i+1;
    
end;

fclose(fid);

[b, ia, ib] = unique(upper(S(:,1:3)),'rows');%Find the unique stations



[c1,i] = setdiff(upper(new_S),upper(S(:,1:8)),'rows');

[c2,i] = setdiff(upper(S(:,1:8)),upper(new_S),'rows');



save D:\Manoj\imag19922004\IMAGfromDVD\missigfiles.mat c1 c2;



%%

%ftp download of missing data

load  D:\Manoj\imag19922004\IMAGfromDVD\missigfiles.mat c1 c2;



localdirectory = 'D:\Manoj\imag19922004\IMAGfromftp\';



f = ftp('ftp.geolab.nrcan.gc.ca','imagnair','data4W0nj');

binary(f);



problem_donwload = 1;



for i = 1:length(c1),%c1 from prevous cell results
    
    try
        
        if str2double(c1(i,4:5)) <80,
            
            addon = '20';
            
        else
            
            addon = '19';
            
        end;
        
        
        
        remotedirectory = ['/var/ftp/data/definitive/' addon c1(i,4:5) '/' lower(c1(i,1:3)) '/'];
        
        filename = [lower(c1(i,:)) '.zip'];
        
        
        
        
        
        
        
        S = dir([localdirectory filename]);
        
        
        
        if ~isempty(S),
            
            sss = S.bytes ;
            
        else
            
            sss = 300000;
            
        end;
        
        
        
        if isempty(S) || sss < 100,
            
            hh = dir(f,[remotedirectory filename]);
            
            if ~isempty(hh),
                
                
                
                
                
                fprintf('Downloading %s %s %d\n',remotedirectory,filename,i);
                
                try
                    
                    cd (f, remotedirectory);
                    
                    mget(f, filename, localdirectory);
                    
                    hh = dir(f,[remotedirectory filename]);
                    
                    fprintf('Downloaded %s %s File Size = %d\n',remotedirectory,filename,hh.bytes);
                    
                    if (hh.bytes < 10),
                        
                        problem_file_string(problem_download,:) = filename;
                        
                        problem_download = problem_download + 1;
                        
                    end;
                    
                catch me;
                    
                    fprintf('Error at %s %s \n',remotedirectory,filename);
                    
                end;
                
            else
                
                fprintf('Data not available : %s %d \n',[remotedirectory filename],i);
                
                close(f);
                
                f = ftp('ftp.geolab.nrcan.gc.ca','imagnair','data4W0nj');
                
                
                
            end;
            
        end;
        
    catch me1;
        
        fprintf('Error somwhere, retrying : %s %d \n',[remotedirectory filename],i);
        
        
        
        
        
    end;
    
end;



close(f);





%%

% Just copy all the files into proper directories



source_dir = 'D:\Manoj\imag19922004\IMAGfromftp\extra\';

copy_dir = 'D:\Manoj\imag19922004\IMAGfromDVD\';

S  = dir('D:\Manoj\imag19922004\IMAGfromftp\extra\*.zip');



for i = 1:length(S),
    
    filestring(i,:) = S(i).name;
    
end;



i_count = 1;

for i = 1:length(filestring),
    
    
    
    ss = dir([copy_dir 'MAG' filestring(i,4:5) '\' filestring(i,1:3) '\']);
    
    if isempty(ss),
        
        [status, result] = mkdir([copy_dir 'MAG' filestring(i,4:5) '\'],filestring(i,1:3));
        
    end;
    
    ss = dir([copy_dir 'MAG' filestring(i,4:5) '\' filestring(i,1:3) '\' filestring(i,:)]);
    
    
    
    
    
    if isempty(ss),
        
        %display([copy_dir 'MAG' filestring(i,4:5) '\' filestring(i,1:3) '\' filestring(i,:)]);
        
        [status, result] = system(['cp ' source_dir filestring(i,:) ' ' copy_dir 'MAG' filestring(i,4:5) '\' filestring(i,1:3) '\' filestring(i,:) ]);
        
    else
        
        display('OK');
        
        i_count = i_count + 1;
        
    end;
    
    
    
end;



%% download data from WDC edin ftp site



load  D:\Manoj\imag19922004\IMAGfromDVD\missigfiles.mat c1 c2; %c1 - stations that are available at wdc and not avi at intermag



[b, ia, ib] = unique(upper(c1(:,1:5)),'rows');%Find the unique station years



%ftp://ftp.nmh.ac.uk/wdc/obsdata/1minval/2004/



localdirectory = 'D:\Manoj\imag19922004\WDCfromftp\';



f = ftp('ftp.nmh.ac.uk');

binary(f);



problem_donwload = 1;



for i = 1:length(b),%c1 from prevous cell results
    
    try
        
        if str2double(b(i,4:5)) <80,
            
            addon = '20';
            
        else
            
            addon = '19';
            
        end;
        
        
        
        remotedirectory = ['/pub/wdc/obsdata/1minval/' addon b(i,4:5) '/'] ;
        
        filename = [lower(b(i,1:3)) addon b(i,4:5) 'm.zip'];
        
        
        
        
        
        S = dir([localdirectory filename]);
        
        
        
        if ~isempty(S),
            
            sss = S.bytes ;
            
        else
            
            sss = 300000;
            
        end;
        
        
        
        if isempty(S) || sss < 100,
            
            hh = dir(f,[remotedirectory filename]);
            
            if ~isempty(hh),
                
                
                
                
                
                fprintf('Downloading %s %s %d\n',remotedirectory,filename,i);
                
                try
                    
                    cd (f, remotedirectory);
                    
                    mget(f, filename, localdirectory);
                    
                    hh = dir(f,[remotedirectory filename]);
                    
                    fprintf('Downloaded %s %s File Size = %d\n',remotedirectory,filename,hh.bytes);
                    
                    if (hh.bytes < 10),
                        
                        problem_file_string(problem_download,:) = filename;
                        
                        problem_download = problem_download + 1;
                        
                    end;
                    
                catch me;
                    
                    fprintf('Error at %s %s \n',remotedirectory,filename);
                    
                end;
                
            else
                
                fprintf('Data not available : %s %d \n',[remotedirectory filename],i);
                
                close(f);
                
                f = ftp('ftp.nmh.ac.uk');
                
                
                
            end;
            
        end;
        
    catch me1;
        
        fprintf('Error somwhere, retrying : %s %d \n',[remotedirectory filename],i);
        
        
        
        
        
    end;
    
end;



close(f);







%% unzip, cat and read the 1 min wdc yearly ascii (iag 2002) file



localdirectory = 'D:\Manoj\imag19922004\WDCfromftp\new\';



unzipdirectory = 'D:\Manoj\imag19922004\WDCfromftp\new\';



S = dir([localdirectory '*.zip']);



for i = 1:length(S),
    
    unzip([localdirectory S(i).name],unzipdirectory);
    
    SS = dir([unzipdirectory S(i).name([1,2,3,6,7])  '*.wdc']);
    
    if ~isempty(SS),
        
        SSS = dir([unzipdirectory S(i).name([1,2,3,6,7])  'all.wdc' ]);
        
        if isempty(SSS),
            
            [status, result] = system(['copy ' unzipdirectory S(i).name([1,2,3,6,7])  '*.wdc ' unzipdirectory S(i).name([1,2,3,6,7])  'all.wdc' ]);
            
        end;
        
    end;
    
end;



%% Add wdc edin 1 min data to the netcdf archive





%%

% reading 1 min binaru imag data



S  = dir('D:\Manoj\imag19922004\WDCfromftp\matlabdatafiles\*.mat');

localdirectory = 'D:\Manoj\imag19922004\WDCfromftp\matlabdatafiles\';

start_date = datenum(1995,1,1);

a = sum(daysinyear([1995:2006]))*24*60;



for i = 1:length(S),
    
    filestring(i,:) = S(i).name;
    
end;



% for i= 1:length(S),

%     filestring(i,:) = cell2str(S(i));

% end;



[b, ia, ib] = unique(filestring(:,1:3),'rows');%Find the unique stations



for i = 1:length(b),
    
    
    
    nc_fname = ['C:\Manoj\projects\gary\minute_data\' b(i,:) '.nc'];
    
    S1 = dir(nc_fname);
    
    
    
    if isempty(S1),
        
        ncid = netcdf.create(nc_fname,'NC_NOCLOBBER');
        
        data_time_axis_dim = netcdf.defDim(ncid,'data_time_axis_dim',a);
        
        X_ID = netcdf.defVar(ncid,'X','int', data_time_axis_dim);
        
        Y_ID = netcdf.defVar(ncid,'Y','int', data_time_axis_dim);
        
        Z_ID = netcdf.defVar(ncid,'Z','int', data_time_axis_dim);
        
        netcdf.putAtt(ncid,X_ID,'_FillValue',int32(999999));
        
        netcdf.putAtt(ncid,Y_ID,'_FillValue',int32(999999));
        
        netcdf.putAtt(ncid,Z_ID,'_FillValue',int32(999999));
        
        netcdf.endDef(ncid);
        
    else
        
        ncid = netcdf.open(nc_fname,'NC_WRITE');
        
        X_ID = netcdf.inqVarID(ncid,'X');
        
        Y_ID = netcdf.inqVarID(ncid,'Y');
        
        Z_ID = netcdf.inqVarID(ncid,'Z');
        
    end;
    
    
    
    L = ib == i; %select a station
    
    filesubset = filestring(L,:);
    
    
    
    
    
    for j = 1:sum(L),
        
        
        
        
        
        matfname = [localdirectory filesubset(j,:)];
        
        
        
        S = dir(matfname);
        
        
        
        if isempty(S),
            
            %no wdc data found
            
            display('error');
            
            break
            
        end;
        
        
        
        eval(['load ' matfname]);
        
        
        
        if ~isempty(X_fday) && max(abs(diff(X_fday))) < 1.1/(24*60),
            
            
            
            X_ts(isnan(X_ts)) = 99999.9;    %no data (this will get converted to 999999 by int32(99999.9*10)
            
            fprintf('file name: %s\n',matfname);
            
            where_to_store = ( (X_fday(1)-(30/(24*3600))) - start_date ) * 1440 ; %30/24*3600 is since X-fday is centered on 30 sec (middle of a minite)
            
            if round(where_to_store) ~= where_to_store,
                
                where_to_store = round(where_to_store);
                
                display(['Warning. Rounded the start time stamp']);
                
            end;
            
            length_of_record = length(X_fday);
            
            
            
            %netcdf.putVar(ncid,varid,start,count,data)
            
            %multiplied by 10 to preserve the decimal points
            
            netcdf.putVar(ncid, X_ID, where_to_store, length_of_record,  int32 (X_ts .* 10) );
            
        else
            
            display(['Error X ' matfname]);
            
            
            
        end;
        
        
        
        if ~isempty(X_fday) && max(abs(diff(X_fday))) < 1.1/(24*60),
            
            
            
            Y_ts(isnan(Y_ts)) = 99999.9;    %no data (this will get converted to 999999 by int32(99999.9*10)
            
            fprintf('file name: %s\n',matfname);
            
            where_to_store = ( (X_fday(1)-(30/(24*3600))) - start_date ) * 1440 ; %30/24*3600 is since X-fday is centered on 30 sec (middle of a minite)
            
            if round(where_to_store) ~= where_to_store,
                
                where_to_store = round(where_to_store);
                
                display(['Warning. Rounded the start time stamp']);
                
            end;
            
            length_of_record = length(X_fday);
            
            
            
            %netcdf.putVar(ncid,varid,start,count,data)
            
            %multiplied by 10 to preserve the decimal points
            
            netcdf.putVar(ncid, Y_ID, where_to_store, length_of_record,  int32 (Y_ts .* 10) );
            
        else
            
            display(['Error Y ' matfname]);
            
            
            
        end;
        
        
        
        if ~isempty(X_fday) && max(abs(diff(X_fday))) < 1.1/(24*60),
            
            
            
            Z_ts(isnan(Z_ts)) = 99999.9;    %no data (this will get converted to 999999 by int32(99999.9*10)
            
            fprintf('file name: %s\n',matfname);
            
            where_to_store = ( (X_fday(1)-(30/(24*3600))) - start_date ) * 1440 ; %30/24*3600 is since X-fday is centered on 30 sec (middle of a minite)
            
            if round(where_to_store) ~= where_to_store,
                
                where_to_store = round(where_to_store);
                
                display(['Warning. Rounded the start time stamp']);
                
            end;
            
            length_of_record = length(X_fday);
            
            
            
            %netcdf.putVar(ncid,varid,start,count,data)
            
            %multiplied by 10 to preserve the decimal points
            
            netcdf.putVar(ncid, Z_ID, where_to_store, length_of_record,  int32 (Z_ts .* 10) );
            
        else
            
            display(['Error Z ' matfname]);
            
            
            
        end;
        
        
        
        
        
        
        
        
        
        
        
    end;
    
    
    
    netcdf.close(ncid);
    
    
    
end;











%% Add wdc edin 1 min data to the netcdf archive





%%

a = sum(daysinyear([1995:2006]))*24*60;

% reading 1 min binaru imag data



S  = dir('D:\Manoj\imag19922004\WDCfromftp\new\*.wdc');

localdirectory = 'D:\Manoj\imag19922004\WDCfromftp\matlabdatafiles\';



for i = 1:length(S),
    
    filestring(i,:) = S(i).name;
    
end;



% for i= 1:length(S),

%     filestring(i,:) = cell2str(S(i));

% end;



[b, ia, ib] = unique(filestring(:,1:3),'rows');%Find the unique stations



for i = 1:length(b),
    
    
    
    nc_fname = ['C:\Manoj\projects\gary\minute_data\' b(i,:) '.nc'];
    
    S1 = dir(nc_fname);
    
    
    
    if ~isempty(S1),
        
        
        
        ncid = netcdf.open(nc_fname,'NC_WRITE');
        
        X_ID = netcdf.inqVarID(ncid,'X');
        
        Y_ID = netcdf.inqVarID(ncid,'Y');
        
        Z_ID = netcdf.inqVarID(ncid,'Z');
        
        
        
        L = ib == i; %select a station
        
        filesubset = filestring(L,:);
        
        
        
        x_data = netcdf.getVar(ncid, X_ID, 0, a);
        
        
        
        x_data = double(x_data)/10;
        
        x_data(x_data==99999.9) = NaN;
        
        plot(datenum(1995,1,1,0,0,30):(1/(24*60)):datenum(2006,12,31,23,59,30),x_data);
        
        hold on;
        
        
        
        for j = 1:sum(L),
            
            
            
            
            
            wdcfname = ['D:\Manoj\imag19922004\WDCfromftp\new\' filesubset(j,1:9) 'wdc'];
            
            
            
            S = dir(wdcfname);
            
            
            
            if isempty(S),
                
                %no wdc data found
                
                display('error');
                
                break
                
            end;
            
            
            
            [X_fday,Y_fday,Z_fday,X_ts,Y_ts,Z_ts] = obs_1_min_read(wdcfname);
            
            
            
            plot(X_fday,X_ts,'r');
            
            hold on;
            
            display(wdcfname);
            
            %eval(['save ' localdirectory filesubset(j,1:9) '.mat' ' X_ts Y_ts Z_ts X_fday']);
            
            
            
        end;
        
        
        
        title(['Minute data 1995 - 2006 ' nc_fname(end-5:end) ' ' b(i,:) ]);
        
        datetick('x','yy');
        
        axis([datenum(1995,1,1) datenum(2006,12,31) -inf inf]);
        
        
        
        netcdf.close(ncid);
        
        %pause;
        
        saveas(gcf,[localdirectory filesubset(j,1:3) '.pdf'],'pdf');
        
        clf;
        
    end;
    
    
    
end;



%% To read the WDC files stored as MATLAB data, compare with the

% data in netcdf format



a = sum(daysinyear([1995:2006]))*24*60;

% reading 1 min binaru imag data



S  = dir('D:\Manoj\imag19922004\WDCfromftp\unzippedfiles\*.wdc');

localdirectory = 'D:\Manoj\imag19922004\WDCfromftp\matlabdatafiles\';



for i = 1:length(S),
    
    filestring(i,:) = S(i).name;
    
end;



% for i= 1:length(S),

%     filestring(i,:) = cell2str(S(i));

% end;



[b, ia, ib] = unique(filestring(:,1:3),'rows');%Find the unique stations



for i = 1:length(b),
    
    
    
    nc_fname = ['C:\Manoj\projects\gary\minute_data\' b(i,:) '.nc'];
    
    S1 = dir(nc_fname);
    
    
    
    if ~isempty(S1),
        
        
        
        ncid = netcdf.open(nc_fname,'NC_WRITE');
        
        X_ID = netcdf.inqVarID(ncid,'X');
        
        Y_ID = netcdf.inqVarID(ncid,'Y');
        
        Z_ID = netcdf.inqVarID(ncid,'Z');
        
        
        
        L = ib == i; %select a station
        
        filesubset = filestring(L,:);
        
        
        
        z_data = netcdf.getVar(ncid, Y_ID, 0, a);
        
        
        
        z_data = double(z_data)/10;
        
        z_data(z_data==99999.9) = NaN;
        
        plot(datenum(1995,1,1,0,0,30):(1/(24*60)):datenum(2006,12,31,23,59,30),z_data);
        
        title(nc_fname(end-7:end));
        
        hold on;
        
        
        
        for j = 1:sum(L),
            
            
            
            
            
            wdcfname = [localdirectory filesubset(j,1:9) '.mat'];
            
            
            
            S = dir(wdcfname);
            
            
            
            if isempty(S),
                
                %no wdc data found
                
                display('error');
                
                break
                
            end;
            
            
            
            eval(['load ' wdcfname]);
            
            
            
            plot(X_fday,Z_ts,'r');
            
            hold on;
            
            display(wdcfname);
            
            
            
        end;
        
        
        
        datetick
        
        netcdf.close(ncid);
        
        %pause;
        
        saveas(gcf,[localdirectory filesubset(j,1:3) '_z.pdf'],'pdf');
        
        clf;
        
    end;
    
    
    
end;





%% Correct the data





nc_fname = 'C:\Manoj\projects\gary\minute_data\ABG.nc';

a = sum(daysinyear([1995:2006]))*24*60;

start_date = datenum(1995,1,1);



S1 = dir(nc_fname);



if ~isempty(S1),
    
    ncid = netcdf.open(nc_fname,'NC_WRITE');
    
    X_ID = netcdf.inqVarID(ncid,'X');
    
    Y_ID = netcdf.inqVarID(ncid,'Y');
    
    Z_ID = netcdf.inqVarID(ncid,'Z');
    
end;



x_data = netcdf.getVar(ncid, X_ID, 0, a);

x_data = double(x_data)/10;

x_data(x_data==99999.9) = NaN;



y_data = netcdf.getVar(ncid, Y_ID, 0, a);

y_data = double(y_data)/10;

y_data(y_data==99999.9) = NaN;



z_data = netcdf.getVar(ncid, Z_ID, 0, a);

z_data = double(z_data)/10;

z_data(z_data==99999.9) = NaN;



%VSS

% L = x_data > 1.9e4;

%

% x_data(L) = NaN;

% y_data(L) = NaN;

% z_data(L) = NaN;

%

% L = abs(z_data) < 1.25e4;

%

% x_data(L) = NaN;

% y_data(L) = NaN;

% z_data(L) = NaN;





%TUC

% a1 = sum(daysinyear([1995]))*24*60;%break here

% z_data(1:a1) = z_data(1:a1)+diff(z_data(a1:a1+1));

% y_data(1:a1) = y_data(1:a1)+diff(y_data(a1:a1+1));

%z_data(662372) = NaN; %spike



%MBC

%ss=find(x_data==max(x_data));

%x_data(ss) = NaN;%wrong no data value,. removed



%LMM

% L = abs(z_data) > 2.6e4;

% z_data(L) = NaN;

% x_data(L) = NaN;

% y_data(L) = NaN;

%L = y_data > -4000;

% z_data(L) = NaN;

% x_data(L) = NaN;

% y_data(L) = NaN;



%HUA

% L = x_data > 5e4;

% z_data(L) = NaN;

% x_data(L) = NaN;

% y_data(L) = NaN;



%GUA

% L = z_data > 8000;

% sum(L)

% z_data(L)

% z_data(L) = NaN;

% x_data(L) = NaN;

% y_data(L) = NaN;

%



%BJN

% L = x_data > 1.5e4;

% z_data(L) = NaN;

% x_data(L) = NaN;

% y_data(L) = NaN;





%ABG

% L = x_data > 3.9e4;

% z_data(L) = NaN;

% x_data(L) = NaN;

% y_data(L) = NaN;

% L = y_data > 2e3;

% z_data(L) = NaN;

% x_data(L) = NaN;

% y_data(L) = NaN;





y_data(isnan(y_data)) = 99999.9;    %no data (this will get converted to 999999 by int32(99999.9*10)

where_to_store = 0;

length_of_record = a;

netcdf.putVar(ncid, Y_ID, where_to_store, length_of_record,  int32 (y_data .* 10) );



x_data(isnan(x_data)) = 99999.9;    %no data (this will get converted to 999999 by int32(99999.9*10)

where_to_store = 0;

length_of_record = a;

netcdf.putVar(ncid, X_ID, where_to_store, length_of_record,  int32 (x_data .* 10) );



z_data(isnan(z_data)) = 99999.9;    %no data (this will get converted to 999999 by int32(99999.9*10)

where_to_store = 0;

length_of_record = a;

netcdf.putVar(ncid, Z_ID, where_to_store, length_of_record,  int32 (z_data .* 10) );



netcdf.close(ncid);





%% update of 2012



%%



datadirectory = '/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/work/intermagnet/ftp.intermagnet.org/definitive/';



fid = fopen('/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/work/1_min_log_error.txt','wt');



S = dir(datadirectory);

S(1:2) = [];

start_date = datenum(2008,1,1);

a = sum(daysinyear(2008:2012))*24*60;

for i= 1:length(S),
    
    filestring(i,:) = S(i).name;
    
end;



[b, ia, ib] = unique(filestring(:,1:3),'rows');%Find the unique stations



for i = 1:length(b),
    
    
    
    nc_fname = ['/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/work/nc/' b(i,:) '.nc'];
    
    
    
    S1 = dir(nc_fname);
    
    
    
    if isempty(S1),
        
        ncid = netcdf.create(nc_fname,'NC_NOCLOBBER');
        
        data_time_axis_dim = netcdf.defDim(ncid,'data_time_axis_dim',a);
        
        X_ID = netcdf.defVar(ncid,'X','int', data_time_axis_dim);
        
        Y_ID = netcdf.defVar(ncid,'Y','int', data_time_axis_dim);
        
        Z_ID = netcdf.defVar(ncid,'Z','int', data_time_axis_dim);
        
        netcdf.putAtt(ncid,X_ID,'_FillValue',int32(999999));
        
        netcdf.putAtt(ncid,Y_ID,'_FillValue',int32(999999));
        
        netcdf.putAtt(ncid,Z_ID,'_FillValue',int32(999999));
        
        netcdf.endDef(ncid);
        
    else
        
        ncid = netcdf.open(nc_fname,'NC_WRITE');
        
        X_ID = netcdf.inqVarID(ncid,'X');
        
        Y_ID = netcdf.inqVarID(ncid,'Y');
        
        Z_ID = netcdf.inqVarID(ncid,'Z');
        
    end;
    
    
    
    L = ib == i; %select a station
    
    filesubset = filestring(L,:);
    
    
    
    
    
    for j = 1:sum(L),
        
        
        
        
        
        binfname = [datadirectory filesubset(j,1:8) '.bin'];
        
        S = dir(binfname);
        
        
        
        if isempty(S),
            
            zipfilename = [datadirectory filesubset(j,1:8) '.zip'];
            
            unzip(zipfilename, datadirectory);
            
        end;
        
        
        
        
        
        
        
        data_struct = INTERMAG(binfname, 'XYZF');
        
        
        
        if ~isempty(data_struct) && (max(abs(diff(data_struct.MINUTELY(:,1)))) < 1.1/(24*60)),
            
            
            
            data_struct.MINUTELY(isnan(data_struct.MINUTELY)) = 99999.9;    %no data (this will get converted to 999999 by int32(99999.9*10)
            
            fprintf('file name: %s\n',filesubset(j,1:8));
            
            where_to_store = ( data_struct.MINUTELY(1,1) - start_date ) * 1440 ;
            
            length_of_record = length(data_struct.MINUTELY(:,1));
            
            
            
            %netcdf.putVar(ncid,varid,start,count,data)
            
            %multiplied by 10 to preserve the decimal points
            
            netcdf.putVar(ncid, X_ID, where_to_store, length_of_record,  int32 (data_struct.MINUTELY(:,2) .* 10) );
            
            netcdf.putVar(ncid, Y_ID, where_to_store, length_of_record,  int32 (data_struct.MINUTELY(:,3) .* 10) );
            
            netcdf.putVar(ncid, Z_ID, where_to_store, length_of_record,  int32 (data_struct.MINUTELY(:,4) .* 10) );
            
            fprintf(fid,'%s             OK \n',filesubset(j,1:8));
            
        else
            
            fprintf(fid,'%s             ERROR \n',filesubset(j,1:8));
            
            
            
        end;
        
        
        
    end;
    
    
    
    netcdf.close(ncid);
    
    
    
end;

fclose(fid);



%% Now add the preliminary data, if needed



datadirectory = '/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/work/intermagnet/ftp.intermagnet.org/preliminary/';
S = dir('/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/work/intermagnet/ftp.intermagnet.org/preliminary/*.min');
fid = fopen('/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/work/1_min_intermagnet_prelim_log_error.txt','at');





% Use this command find ./ -type f -exec cp '{}' ./ \; to copy all the

% files to one directory



more off;





for i= 1:length(S),
    
    filestring(i,:) = S(i).name;
    
end;



usgsobs = ['brw';'bou';'cmo';'ded';'frd';'frn';'gua';'hon';'new';'sjg';'shu';'sit';'bsl';'tuc'];

[b, ia, ib] = unique(filestring(:,1:3),'rows');%Find the unique stations

%select only USGS observatories

[b_usgs, i_aa, i_bb] = intersect(b,usgsobs,'rows');



a = (datenum(2012,12,31) - datenum(2010,1,1) + 1) * 24 * 60;

start_date = datenum(2010,1,1);

%%

%for i = 1:length(b_usgs),

for kk = 1:length(i_aa), % only usgs observatories
    
    i = i_aa(kk);
    
    nc_fname = ['/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/work/intermagnet/ftp.intermagnet.org/preliminary_2010_2012_nc/' b(i,:) '_prelim.nc'];
    
    
    
    S1 = dir(nc_fname);
    
    
    
    if isempty(S1),
        
        ncid = netcdf.create(nc_fname,'NC_NOCLOBBER');
        
        data_time_axis_dim = netcdf.defDim(ncid,'data_time_axis_dim',a);
        
        X_ID = netcdf.defVar(ncid,'X','int', data_time_axis_dim);
        
        Y_ID = netcdf.defVar(ncid,'Y','int', data_time_axis_dim);
        
        Z_ID = netcdf.defVar(ncid,'Z','int', data_time_axis_dim);
        
        netcdf.putAtt(ncid,X_ID,'_FillValue',int32(999999));
        
        netcdf.putAtt(ncid,Y_ID,'_FillValue',int32(999999));
        
        netcdf.putAtt(ncid,Z_ID,'_FillValue',int32(999999));
        
        netcdf.endDef(ncid);
        
    else
        
        ncid = netcdf.open(nc_fname,'NC_WRITE');
        
        X_ID = netcdf.inqVarID(ncid,'X');
        
        Y_ID = netcdf.inqVarID(ncid,'Y');
        
        Z_ID = netcdf.inqVarID(ncid,'Z');
        
    end;
    
    
    
    L = ib == i; %select a station
    
    filesubset = filestring(L,:);
    
    
    
    
    
    for j = 1:sum(L),
        
        
        
        
        
        minfname = [datadirectory filesubset(j,:)];
        
        
        
        
        
        data_struct = intermagnet(minfname, 'XYZF');
        
        
        
        if (~isempty(data_struct))
            
            data_struct.DATA(isnan(data_struct.DATA)) = 99999.9;    %no data (this will get converted to 999999 by int32(99999.9*10)
            
            if(strcmp(data_struct.CORD,'XYZF')),
                
                display('HDZF to XYZF conversion');
                
            end;
            
            fprintf('file name: %s\n',filesubset(j,:));
            
            where_to_store = ( datenum(filesubset(j,4:11),'yyyymmdd') - start_date ) * 1440 ;
            
            %netcdf.putVar(ncid,varid,start,count,data)
            
            % multiplied by 10 to preserve the decimal points
            
            netcdf.putVar(ncid, X_ID, where_to_store, 1440, int32( data_struct.DATA(:,1) .* 10) );
            
            netcdf.putVar(ncid, Y_ID, where_to_store, 1440, int32( data_struct.DATA(:,2) .* 10) );
            
            netcdf.putVar(ncid, Z_ID, where_to_store, 1440, int32( data_struct.DATA(:,3) .* 10) );
            
            
            
        else,
            
            fprintf(fid,'%s             ERROR \n',filesubset(j,:));
            
        end;
        
        
        
    end;
    
    
    
    netcdf.close(ncid);
    
    
    
    
    
end;



%% Add wdc edin 1 min data to the netcdf archive





%%



datadirectory = '/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/work/WDCEdin/';



fid = fopen('/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/work/1_min_wdc_edin_log_error.txt','wt');



S = dir('/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/work/WDCEdin/*.wdc');



start_date = datenum(2005,1,1); % 2005 since LRV

%a = sum(daysinyear(2008:2011))*24*60;

a = (datenum(2011,12,31) - datenum(2005,1,1) + 1) * 24 * 60;



for i = 1:length(S),
    
    filestring(i,:) = S(i).name;
    
end;





[b, ia, ib] = unique(filestring(:,1:3),'rows');%Find the unique stations



for i = 14:length(b),
    
    
    
    nc_fname = ['/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/work/WDCEdin/nc/' b(i,:) '.nc'];
    
    S1 = dir(nc_fname);
    
    
    
    if isempty(S1),
        
        ncid = netcdf.create(nc_fname,'NC_NOCLOBBER');
        
        data_time_axis_dim = netcdf.defDim(ncid,'data_time_axis_dim',a);
        
        X_ID = netcdf.defVar(ncid,'X','int', data_time_axis_dim);
        
        Y_ID = netcdf.defVar(ncid,'Y','int', data_time_axis_dim);
        
        Z_ID = netcdf.defVar(ncid,'Z','int', data_time_axis_dim);
        
        netcdf.putAtt(ncid,X_ID,'_FillValue',int32(999999));
        
        netcdf.putAtt(ncid,Y_ID,'_FillValue',int32(999999));
        
        netcdf.putAtt(ncid,Z_ID,'_FillValue',int32(999999));
        
        netcdf.endDef(ncid);
        
    else
        
        ncid = netcdf.open(nc_fname,'NC_WRITE');
        
        X_ID = netcdf.inqVarID(ncid,'X');
        
        Y_ID = netcdf.inqVarID(ncid,'Y');
        
        Z_ID = netcdf.inqVarID(ncid,'Z');
        
    end;
    
    
    
    L = ib == i; %select a station
    
    filesubset = filestring(L,:);
    
    
    
    for j = 1:sum(L),
        
        
        
        
        
        wdcfname = [datadirectory filesubset(j,:)];
        
        
        
        S2 = dir(wdcfname);
        
        
        
        if isempty(S2),
            
            %no wdc data found
            
            display('error');
            
            break
            
        end;
        
        
        
        [X_fday,Y_fday,Z_fday,X_ts,Y_ts,Z_ts] = obs_1_min_read(wdcfname);
        
        
        
        % plot(X_fday,X_ts,'r');
        
        % hold on;
        
        % display(wdcfname);
        
        %eval(['save ' localdirectory filesubset(j,1:9) '.mat' ' X_ts Y_ts Z_ts X_fday']);
        
        
        
        if (length(X_fday) == length(Y_fday) &&  length(X_fday)  == length(Z_fday) && ~any(diff(X_fday) > 1.05/1440))
            
            
            
            
            
            fprintf('file name: %s\n',filesubset(j,:));
            
            %where_to_store = ( datenum(filesubset(j,4:11),'yyyymmdd') - start_date ) * 1440 ;
            
            where_to_store = ((X_fday(1) - (30/(3600*24))) - start_date )*1440; % remove 30 seconds
            
            %netcdf.putVar(ncid,varid,start,count,data)
            
            % multiplied by 10 to preserve the decimal points
            
            netcdf.putVar(ncid, X_ID, where_to_store, length(X_ts), int32( X_ts .* 10) );
            
            netcdf.putVar(ncid, Y_ID, where_to_store, length(Y_ts), int32( Y_ts .* 10) );
            
            netcdf.putVar(ncid, Z_ID, where_to_store, length(Z_ts), int32( Z_ts .* 10) );
            
            
            
            subplot(311);
            
            plot(X_fday,X_ts);
            
            title(filesubset(j,:));
            
            hold on;
            
            subplot(312);
            
            plot(Y_fday,Y_ts);hold on;
            
            subplot(313);
            
            plot(Z_fday,Z_ts);hold on;
            
            
            
        else,
            
            
            
            fprintf(fid,'ERROR  %s\n', filesubset(j,:));
            
            
            
        end;
        
        
        
    end;
    
    
    
    netcdf.close(ncid);
    
    datetick('x','yy/mm','keeplimits');
    
    saveas(gcf,['/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/work/WDCEdin/plots/' b(i,:) '_wdc'],'png');
    
    close all;
    
    
    
end;





%% Now add the WDC Kyoto data



datadirectory = '/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/work/WDCKyoto/';

S = dir('/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/work/WDCKyoto/*.dat');

fid = fopen('/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/work/1_min_Kyoto_intermagnet_prelim_log_error.txt','at');





% Use this command find ./ -type f -exec cp '{}' ./ \; to copy all the

% files to one directory



more off;





for i= 1:length(S),
    
    filestring(i,:) = S(i).name;
    
end;



[b, ia, ib] = unique(filestring(:,1:3),'rows');%Find the unique stations



a = (datenum(2012,12,31) - datenum(2010,1,1) + 1) * 24 * 60;

start_date = datenum(2010,1,1);



for i = 1:length(b),
    
    
    
    nc_fname = ['/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/work/WDCKyoto/nc/' b(i,:) '_prelim.nc'];
    
    
    
    S1 = dir(nc_fname);
    
    
    
    if isempty(S1),
        
        ncid = netcdf.create(nc_fname,'NC_NOCLOBBER');
        
        data_time_axis_dim = netcdf.defDim(ncid,'data_time_axis_dim',a);
        
        X_ID = netcdf.defVar(ncid,'X','int', data_time_axis_dim);
        
        Y_ID = netcdf.defVar(ncid,'Y','int', data_time_axis_dim);
        
        Z_ID = netcdf.defVar(ncid,'Z','int', data_time_axis_dim);
        
        netcdf.putAtt(ncid,X_ID,'_FillValue',int32(999999));
        
        netcdf.putAtt(ncid,Y_ID,'_FillValue',int32(999999));
        
        netcdf.putAtt(ncid,Z_ID,'_FillValue',int32(999999));
        
        netcdf.endDef(ncid);
        
    else
        
        ncid = netcdf.open(nc_fname,'NC_WRITE');
        
        X_ID = netcdf.inqVarID(ncid,'X');
        
        Y_ID = netcdf.inqVarID(ncid,'Y');
        
        Z_ID = netcdf.inqVarID(ncid,'Z');
        
    end;
    
    
    
    L = ib == i; %select a station
    
    filesubset = filestring(L,:);
    
    
    
    
    
    for j = 1:sum(L),
        
        
        
        
        
        minfname = [datadirectory filesubset(j,:)];
        
        
        
        
        
        data_struct = intermagnet(minfname, 'XYZF');
        
        
        
        if (~isempty(data_struct))
            
            data_struct.DATA(isnan(data_struct.DATA)) = 99999.9;    %no data (this will get converted to 999999 by int32(99999.9*10)
            
            if(strcmp(data_struct.CORD,'XYZF')),
                
                display('HDZF to XYZF conversion');
                
            end;
            
            fprintf('file name: %s\n',filesubset(j,:));
            
            where_to_store = ( datenum(filesubset(j,4:11),'yyyymmdd') - start_date ) * 1440 ;
            
            %netcdf.putVar(ncid,varid,start,count,data)
            
            % multiplied by 10 to preserve the decimal points
            
            netcdf.putVar(ncid, X_ID, where_to_store, 1440, int32( data_struct.DATA(:,1) .* 10) );
            
            netcdf.putVar(ncid, Y_ID, where_to_store, 1440, int32( data_struct.DATA(:,2) .* 10) );
            
            netcdf.putVar(ncid, Z_ID, where_to_store, 1440, int32( data_struct.DATA(:,3) .* 10) );
            
            
            
        else,
            
            fprintf(fid,'%s             ERROR \n',filesubset(j,:));
            
        end;
        
        
        
    end;
    
    
    
    netcdf.close(ncid);
    
    
    
    
    
end;



%% Update of 2012 WDC Kyoto ADD KNZ





start_date = datenum(2008,1,1);

a = (datenum(2012,12,31) - datenum(2008,1,1) + 1) * 24 * 60;

time_array_min = (datenum(2008,1,1,0,0,30): (1/(24*60)): datenum(2012,12,31,23,59,30))';







nc_fname = ['/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/work/WDCKyoto/nc/KNZ_prelim.nc'];



S1 = dir(nc_fname);



if isempty(S1),
    
    ncid = netcdf.create(nc_fname,'NC_NOCLOBBER');
    
    data_time_axis_dim = netcdf.defDim(ncid,'data_time_axis_dim',a);
    
    X_ID = netcdf.defVar(ncid,'X','int', data_time_axis_dim);
    
    Y_ID = netcdf.defVar(ncid,'Y','int', data_time_axis_dim);
    
    Z_ID = netcdf.defVar(ncid,'Z','int', data_time_axis_dim);
    
    netcdf.putAtt(ncid,X_ID,'_FillValue',int32(999999));
    
    netcdf.putAtt(ncid,Y_ID,'_FillValue',int32(999999));
    
    netcdf.putAtt(ncid,Z_ID,'_FillValue',int32(999999));
    
    netcdf.endDef(ncid);
    
else
    
    ncid = netcdf.open(nc_fname,'NC_WRITE');
    
    X_ID = netcdf.inqVarID(ncid,'X');
    
    Y_ID = netcdf.inqVarID(ncid,'Y');
    
    Z_ID = netcdf.inqVarID(ncid,'Z');
    
end;





load /nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/work/WDCKyoto/KNZ_2010.mat



if (~isempty(data_struct))
    
    data_struct.DATA(isnan(data_struct.DATA)) = 99999.9;
    
    
    
    where_to_store = ( data_struct.FDAY(1) - start_date ) * 1440 ;
    
    
    
    netcdf.putVar(ncid, X_ID, where_to_store, length(data_struct.DATA(:,1)), int32( data_struct.DATA(:,1) .* 10) );
    
    netcdf.putVar(ncid, Y_ID, where_to_store, length(data_struct.DATA(:,2)), int32( data_struct.DATA(:,2) .* 10) );
    
    netcdf.putVar(ncid, Z_ID, where_to_store, length(data_struct.DATA(:,3)), int32( data_struct.DATA(:,3) .* 10) );
    
end;



load /nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/work/WDCKyoto/KNZ_2011.mat



if (~isempty(data_struct))
    
    data_struct.DATA(isnan(data_struct.DATA)) = 99999.9;
    
    
    
    where_to_store = ( data_struct.FDAY(1) - start_date ) * 1440 ;
    
    
    
    netcdf.putVar(ncid, X_ID, where_to_store, length(data_struct.DATA(:,1)), int32( data_struct.DATA(:,1) .* 10) );
    
    netcdf.putVar(ncid, Y_ID, where_to_store, length(data_struct.DATA(:,2)), int32( data_struct.DATA(:,2) .* 10) );
    
    netcdf.putVar(ncid, Z_ID, where_to_store, length(data_struct.DATA(:,3)), int32( data_struct.DATA(:,3) .* 10) );
    
end;





netcdf.close(ncid);





%%

%% Update of 2012



% Add intermagnet definitive and preliminary. Look for definitive last

% value in the array and add preliminary after that - if any

% Definitive data



start_date = datenum(2008,1,1);

a = (datenum(2012,12,31) - datenum(2008,1,1) + 1) * 24 * 60;

time_array_min = (datenum(2008,1,1,0,0,30): (1/(24*60)): datenum(2012,12,31,23,59,30))';

basedirectory = '/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/work/intermagnet/ftp.intermagnet.org/definitive_2008_2012_nc/';



a_prelim = (datenum(2012,12,31) - datenum(2010,1,1) + 1) * 24 * 60;

time_array_min_prelim = (datenum(2010,1,1,0,0,30): (1/(24*60)): datenum(2012,12,31,23,59,30))';

basedirectory_prelim = '/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/work/intermagnet/ftp.intermagnet.org/preliminary_2010_2012_nc/';



S1 = dir([basedirectory '*.nc']);

S2 = dir([basedirectory_prelim '*.nc']);



if exist('filestring','var')
    
    clear filestring;
    
end;



if exist('filestring_prelims','var')
    
    clear filestring_prelims;
    
end;



for i= 1:length(S1),
    
    filestring(i,:) = S1(i).name(1:3);
    
end;



for i= 1:length(S2),
    
    filestring_prelims(i,:) = S2(i).name(1:3);
    
end;



spike_stn = ['aqu';'brw';'cmo';'ded';'dmc';'eyr';'fur';'hrb';'hyb';'kiv';'lrm';'lyc';'mbo';'nur';'peg';'sit'];



[stn , index_a] = setdiff(filestring_prelims,filestring,'rows','legacy');


usgsobs = ['brw';'bou';'cmo';'ded';'frd';'frn';'gua';'hon';'new';'sjg';'shu';'sit';'bsl';'tuc'];



for i = 1:length(filestring_prelims), % all definitive data
    
    %for i = 20:20
    
    % process only for USGS stations
    
    if ~isempty(intersect(usgsobs,filestring_prelims(i,:),'rows')),
        
        
        nc_fname_prelims = [basedirectory_prelim S2(i).name];
        
        S4 = dir(nc_fname_prelims);
        
        if ~isempty(S4),
            
            ncid_p = netcdf.open(nc_fname_prelims,'NOWRITE');
            
            X_ID = netcdf.inqVarID(ncid_p,'X');
            
            Y_ID = netcdf.inqVarID(ncid_p,'Y');
            
            Z_ID = netcdf.inqVarID(ncid_p,'Z');
            
        end;
        
        
        
        x_data_prelim = netcdf.getVar(ncid_p, X_ID, 0, a_prelim);
        
        x_data_prelim = double(x_data_prelim)/10;
        
        x_data_prelim(x_data_prelim==99999.9) = NaN;
        
        
        
        y_data_prelim = netcdf.getVar(ncid_p, Y_ID, 0, a_prelim);
        
        y_data_prelim = double(y_data_prelim)/10;
        
        y_data_prelim(y_data_prelim==99999.9) = NaN;
        
        
        
        z_data_prelim = netcdf.getVar(ncid_p, Z_ID, 0, a_prelim);
        
        z_data_prelim = double(z_data_prelim)/10;
        
        z_data_prelim(z_data_prelim==99999.9) = NaN;
        
        netcdf.close(ncid_p);
        
        
        
        
        
        % Write the combined data to netCDF files
        
        
        
        
        
        nf_combined_fname = ['/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/work/combinednc/' S2(i).name(1:3) '_combined.nc']
        
        
        
        
        
        S4 = dir(nf_combined_fname);
        
        
        
        if isempty(S4),
            
            ncid_c = netcdf.create(nf_combined_fname,'NC_NOCLOBBER');
            
            data_time_axis_dim = netcdf.defDim(ncid_c,'data_time_axis_dim',a);
            
            X_ID_c = netcdf.defVar(ncid_c,'Magnetic_Field_X','int', data_time_axis_dim);
            
            Y_ID_c = netcdf.defVar(ncid_c,'Magnetic_Field_Y','int', data_time_axis_dim);
            
            Z_ID_c = netcdf.defVar(ncid_c,'Magnetic_Field_Z','int', data_time_axis_dim);
            
            netcdf.putAtt(ncid_c,X_ID,'_FillValue',int32(999999));
            
            netcdf.putAtt(ncid_c,Y_ID,'_FillValue',int32(999999));
            
            netcdf.putAtt(ncid_c,Z_ID,'_FillValue',int32(999999));
            
            netcdf.endDef(ncid_c);
            
        else
            
            ncid_c = netcdf.open(nf_combined_fname,'NC_WRITE');
            
            X_ID_c = netcdf.inqVarID(ncid_c,'Magnetic_Field_X');
            
            Y_ID_c = netcdf.inqVarID(ncid_c,'Magnetic_Field_Y');
            
            Z_ID_c = netcdf.inqVarID(ncid_c,'Magnetic_Field_Z');
            
        end;
        
        
        
        
        
        % conditional editing specific to stations
        
        
        
        if S2(i).name(1:3) == 'ale'
            
            x_data_prelim(1:24*60*366) = NaN; % 2010 data is crap
            
            y_data_prelim(1:24*60*366) = NaN; % 2010 data is crap
            
            z_data_prelim(1:24*60*366) = NaN; % 2010 data is crap
            
        elseif S2(i).name(1:3) == 'cnh'
            
            LLL = x_data_prelim < 0;
            
            x_data_prelim( LLL) = NaN;
            
            y_data_prelim( LLL) = NaN;
            
            z_data_prelim( LLL) = NaN;
            
            clear LLL;
            
        elseif S2(i).name(1:3) == 'mgd'
            
            x_data_prelim(24*60*549:end) = NaN; % > 2011.5 data is crap
            
            y_data_prelim(24*60*549:end) = NaN; % > 2011.5 data is crap
            
            z_data_prelim(24*60*549:end) = NaN; %  > 2011.5 data is crap
            
        elseif S2(i).name(1:3) == 'hua'
            
            x_data_prelim(24*60*931:end) = NaN; % > 2012.5 data is crap
            
            y_data_prelim(24*60*931:end) = NaN; % > 2012.5 data is crap
            
            z_data_prelim(24*60*931:end) = NaN; %  > 2012.5 data is crap
            
        elseif S2(i).name(1:3) == 'knz'
            
            x_data_prelim(1:end) = NaN; % the 2012 file is dummy
            
            y_data_prelim(1:end) = NaN; % the file is dummy
            
            z_data_prelim(1:end) = NaN; % the file is dummy
            
            
            
        end;
        
        
        
        x_factor = 0;
        
        y_factor = 0;
        
        z_factor = 0;
        
        x_prelim = x_data_prelim;
        
        y_prelim = y_data_prelim;
        
        z_prelim = z_data_prelim;
        
        t_prelim = time_array_min_prelim;
        
        
        
        
        
        
        
        
        
        
        
        [pstn,index_c ] = intersect(filestring,filestring_prelims(i,:),'rows','legacy');
        
        
        
        if (~isempty(pstn))
            
            
            
            nc_fname = [basedirectory S1(index_c).name];
            
            
            
            S3 = dir(nc_fname);
            
            if ~isempty(S3),
                
                ncid = netcdf.open(nc_fname,'WRITE');
                
                X_ID = netcdf.inqVarID(ncid,'X');
                
                Y_ID = netcdf.inqVarID(ncid,'Y');
                
                Z_ID = netcdf.inqVarID(ncid,'Z');
                
            end;
            
            
            
            where_to_store = 0 ;
            
            
            
            x_data = netcdf.getVar(ncid, X_ID, 0, a);
            
            netcdf.putVar(ncid_c, X_ID_c, where_to_store, a,int32( x_data) ); % store the definitive part to a new file
            
            x_data = double(x_data)/10;
            
            x_data(x_data==99999.9) = NaN;
            
            
            
            y_data = netcdf.getVar(ncid, Y_ID, 0, a);
            
            netcdf.putVar(ncid_c, Y_ID_c, where_to_store, a,int32( y_data) );
            
            y_data = double(y_data)/10;
            
            y_data(y_data==99999.9) = NaN;
            
            
            
            z_data = netcdf.getVar(ncid, Z_ID, 0, a);
            
            netcdf.putVar(ncid_c, Z_ID_c, where_to_store, a,int32( z_data) );
            
            z_data = double(z_data)/10;
            
            z_data(z_data==99999.9) = NaN;
            
            
            
            netcdf.close(ncid);
            
            
            
            
            
            
            
            %netcdf.putVar(ncid,varid,start,count,data)
            
            % multiplied by 10 to preserve the decimal points
            
            
            
            
            
            
            
            
            
            subplot(311);
            
            plot(time_array_min,x_data);
            
            hold on;
            
            title(S1(index_c).name);
            
            ylabel('X');
            
            datetick;
            
            subplot(312);
            
            plot(time_array_min,y_data);
            
            hold on;
            
            ylabel('Y');
            
            datetick;
            
            subplot(313);
            
            plot(time_array_min,z_data);
            
            hold on;
            
            datetick;
            
            ylabel('Z');
            
            
            
            
            
            % adjust the jumps in the preliminary data
            
            % find the last value's index
            
            
            
            L = isnan(x_data);
            
            k = find(L==0,1,'last');
            
            
            
            if (~isempty(k))
                
                
                
                kk = find(time_array_min_prelim == time_array_min(k+1));
                
                x_median = nanmedian(x_data(k-1000:k));
                
                y_median = nanmedian(y_data(k-1000:k));
                
                z_median = nanmedian(z_data(k-1000:k));
                
                
                
            else,
                
                
                
                kk = 1;
                
                x_median = 0;
                
                y_median = 0;
                
                z_median = 0;
                
            end;
            
            
            
            
            
            x_prelim = x_data_prelim(kk:end);
            
            y_prelim = y_data_prelim(kk:end);
            
            z_prelim = z_data_prelim(kk:end);
            
            t_prelim = time_array_min_prelim(kk:end);
            
            
            
            % find the median value of the first 1000 points in the def data
            
            
            
            
            
            
            
            
            
            % jump removal
            
            
            
            
            
            if (~isempty(k) && length(x_prelim)> 10000 && sum(~isnan(x_prelim)) > 100 )
                
                
                
                x_prelim_median = nanmedian(x_prelim(1:10000));
                
                if isnan(x_prelim_median),
                    
                    x_prelim_median = nanmedian(x_prelim(1:100000));
                    
                end;
                
                y_prelim_median = nanmedian(y_prelim(1:10000));
                
                if isnan(y_prelim_median),
                    
                    y_prelim_median = nanmedian(x_prelim(1:100000));
                    
                end;
                
                z_prelim_median = nanmedian(z_prelim(1:10000));
                
                if isnan(z_prelim_median),
                    
                    z_prelim_median = nanmedian(x_prelim(1:100000));
                    
                end;
                
                
                
            else
                
                x_prelim_median = NaN;
                
                y_prelim_median = NaN;
                
                z_prelim_median = NaN;
                
                
                
            end;
            
            
            % NO JUMP removal
            
              x_prelim_median = NaN;
                
                y_prelim_median = NaN;
                
                z_prelim_median = NaN;
                
            
            
            if all(~isnan([x_median x_prelim_median]))
                
                x_factor = x_median - x_prelim_median;
                
            else,
                
                x_factor = 0;
                
            end;
            
            
            
            if all(~isnan([y_median y_prelim_median]))
                
                y_factor = y_median - y_prelim_median;
                
            else,
                
                y_factor = 0;
                
            end;
            
            if all(~isnan([z_median z_prelim_median]))
                
                z_factor = z_median - z_prelim_median;
                
            else,
                
                z_factor = 0;
                
            end;
            
            
            
        end;
        
        
        
        
        
        % spike detection and removal
        
        
        
        x_prelim = x_data_prelim(kk:end);
        
        t_prelim = time_array_min_prelim(kk:end);
        
        if (sum(~isnan(x_prelim)) > 100 && any(intersect(spike_stn,filestring_prelims(i,:),'rows')))
            
            
            
            LL = isnan(x_prelim);
            
            x_prelim(LL) = [];
            
            t_prelim(LL) = [];
            
            
            
            
            
            
            
            [a3,b3] = robustfit(t_prelim(1:100:end),x_prelim(1:100:end));
            
            
            
            x_prelim = x_data_prelim(kk:end);
            
            t_prelim = time_array_min_prelim(kk:end);
            
            
            
            resid = (x_prelim) - ( a3(1) + a3(2).*t_prelim);
            
            %LL = abs(resid) >500;
            
            LL = abs(resid) > 10*nanstd(resid);
            
            
            
            
            
            x_prelim(LL) = NaN;
            
        end;
        
        
        
        %Y
        
        
        
        y_prelim = y_data_prelim(kk:end);
        
        t_prelim = time_array_min_prelim(kk:end);
        
        if (sum(~isnan(y_prelim))  > 100 && any(intersect(spike_stn,filestring_prelims(i,:),'rows')))
            
            LL = isnan(y_prelim);
            
            y_prelim(LL) = [];
            
            t_prelim(LL) = [];
            
            
            
            [a3,b3] = robustfit(t_prelim(1:100:end),y_prelim(1:100:end));
            
            
            
            y_prelim = y_data_prelim(kk:end);
            
            t_prelim = time_array_min_prelim(kk:end);
            
            
            
            resid = (y_prelim) - ( a3(1) + a3(2).*t_prelim);
            
            %LL = abs(resid) >500;
            
            LL = abs(resid) > 10*nanstd(resid);
            
            
            
            y_prelim(LL) = NaN;
            
        end;
        
        %Z
        
        
        
        
        
        z_prelim = z_data_prelim(kk:end);
        
        t_prelim = time_array_min_prelim(kk:end);
        
        if (sum(~isnan(z_prelim))  > 100 && any(intersect(spike_stn,filestring_prelims(i,:),'rows')))
            
            LL = isnan(z_prelim);
            
            z_prelim(LL) = [];
            
            t_prelim(LL) = [];
            
            
            
            [a3,b3] = robustfit(t_prelim(1:100:end),z_prelim(1:100:end));
            
            
            
            z_prelim = z_data_prelim(kk:end);
            
            t_prelim = time_array_min_prelim(kk:end);
            
                        
            resid = (z_prelim) - ( a3(1) + a3(2).*t_prelim);
            
            %LL = abs(resid) >500;
            
            LL = abs(resid) > 10*nanstd(resid);
            
            z_prelim(LL) = NaN;
            
            %
            
        end;
        
        
        
        
        
        
        
        subplot(311);
        
        
        
        plot(t_prelim,x_prelim + x_factor,'r');
        
        datetick;
        
        hold on;
        
        title(S2(i).name);
        
        ylabel('X');
        
        subplot(312);
        
        
        
        plot(t_prelim,y_prelim + y_factor,'r');
        
        hold on;
        
        ylabel('Y');
        
        datetick;
        
        subplot(313);
        
        
        
        plot(t_prelim,z_prelim + z_factor,'r');
        
        hold on;
        
        datetick;
        
        ylabel('Z');
        
        
        
        
        
        saveas(gcf,['/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/work/intermagnet/plots/' S2(i).name(1:3) '_with_spike'],'png');
        
        display(S2(i).name(1:3));
        
        
        
        if (length(t_prelim) > 0)
            
            
            
            x_prelim = x_prelim + x_factor;
            
            y_prelim = y_prelim + y_factor;
            
            z_prelim = z_prelim + z_factor;
            
            
            
            
            
            x_prelim(isnan(x_prelim)) = 99999.9;
            
            y_prelim(isnan(y_prelim)) = 99999.9;
            
            z_prelim(isnan(z_prelim)) = 99999.9;
            
            
            
            where_to_store = find(time_array_min == t_prelim(1)) - 1;
            
            
            
            netcdf.putVar(ncid_c, X_ID_c, where_to_store, length(x_prelim),int32( x_prelim  * 10 ) );
            
            netcdf.putVar(ncid_c, Y_ID_c, where_to_store, length(y_prelim),int32( y_prelim  * 10 ) );
            
            netcdf.putVar(ncid_c, Z_ID_c, where_to_store, length(z_prelim),int32( z_prelim  * 10 ) );
            
        end;
        
        
        
        netcdf.close(ncid_c);
        
        
        
        close all;
        
        
    end;
    
    
end;









%% Finally Combine the 2012 and 2013 data





S1 = dir('/nfs/satmag_work/mnair/projects/obs_mag_data/update_2011/absolute/*.nc');



S2 = dir('/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/work/combinednc/*.nc');



problem_stn = ['abg';'ebr';'lvv';'ppt';'thy';'vss'];



start_date = datenum(2008,1,1);

a = (datenum(2012,12,31) - datenum(2008,1,1) + 1) * 24 * 60;

time_array_min = (datenum(2008,1,1,0,0,0): (1/(24*60)): datenum(2012,12,31,23,59,0))';

basedirectory = '/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/work/combinednc/';



a_2011 = (datenum(2010,12,31) - datenum(1995,1,1) + 1) * 24 * 60;

time_array_min_2011 = (datenum(1995,1,1,0,0,0): (1/(24*60)): datenum(2010,12,31,23,59,0))';

basedirectory_2011 = '/nfs/satmag_work/mnair/projects/obs_mag_data/update_2011/absolute/';



% director for netCDF files

localdirectory = '/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/absolute/';



S1 = dir([basedirectory '*.nc']);

S2 = dir([basedirectory_2011 '*.nc']);



if exist('filestring','var')
    
    clear filestring;
    
end;



if exist('filestring_2011','var')
    
    clear filestring_2011;
    
end;



for i= 1:length(S1),
    
    filestring(i,:) = upper(S1(i).name(1:3));
    
end;



for i= 1:length(S2),
    
    filestring_2011(i,:) = S2(i).name(1:3);
    
end;



[master_a,master_b] = xlsread('/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/work/Observatory_Master_Metadata.xlsx');



master_b(1,:) = [];



for i= 1:length(master_b),
    
    master_stn(i,:) = upper(cell2str(master_b(i,1)));
    
end;











% get the total stations



[allstn, index_a, index_b] = union(filestring, filestring_2011,'rows');



for i = 1:length(allstn), % 55 GZH
    
    
    
    
    
    %create new netCDF file for writing the data
    
    
    
    [stn , index_a ,index_b ] = intersect(master_stn,allstn(i,:),'rows');
    
    
    
    metadatacellarray = master_b(index_a ,:);
    
    
    
    lat = master_a(index_a ,1);
    
    
    
    long = master_a(index_a ,2);
    
    
    
    alt = master_a(index_a ,3);
    
    
    
    fprintf('Creating %s ',allstn(i,:));
    
    
    
    [ncid_1, X_ID, Y_ID, Z_ID] = Create_Geomag_netCDF(localdirectory, lat, long, alt, metadatacellarray);
    
    
    
    
    
    %read 2011 update data
    
    
    
    ncfname = [basedirectory_2011 allstn(i,:) '_1995_2010.nc'];
    
    
    
    if ~isempty(dir(ncfname))
        
        
        
        fprintf('Adding 2011 data  \n');
        
        
        
        [x_data, y_data, z_data ] = read_geomag_netcdf(ncfname,0,a_2011);
        
        
        
        subplot(311);
        
        plot(time_array_min_2011, x_data);
        
        axis([time_array_min_2011(1) time_array_min(end) -inf inf]);
        
        
        
        datetick;
        
        ylabel('X');
        
        title(allstn(i,:));
        
        hold on;
        
        subplot(312);
        
        plot(time_array_min_2011, y_data);
        
        axis([time_array_min_2011(1) time_array_min(end) -inf inf]);
        
        
        
        datetick;
        
        ylabel('Y');
        
        hold on;
        
        subplot(313);
        
        plot(time_array_min_2011, z_data);
        
        axis([time_array_min_2011(1) time_array_min(end) -inf inf]);
        
        
        
        datetick;
        
        ylabel('Z');
        
        hold on;
        
        
        
        
        
        
        
        x_data(isnan(x_data)) = 99999.9;    %no data (this will get converted to 999999 by int32(99999.9*10)
        
        z_data(isnan(z_data)) = 99999.9;    %no data (this will get converted to 999999 by int32(99999.9*10)
        
        y_data(isnan(y_data)) = 99999.9;    %no data (this will get converted to 999999 by int32(99999.9*10)
        
        
        
        where_to_store = 0;
        
        length_of_record = a_2011;
        
        
        
        netcdf.putVar(ncid_1, X_ID, where_to_store, length_of_record,  int32 (x_data .* 10) );
        
        netcdf.putVar(ncid_1, Y_ID, where_to_store, length_of_record,  int32 (y_data .* 10) );
        
        netcdf.putVar(ncid_1, Z_ID, where_to_store, length_of_record,  int32 (z_data .* 10) );
        
        
        
    end;
    
    
    
    % read 2012 data
    
    
    
    ncfname = [basedirectory lower(allstn(i,:)) '_combined.nc'];
    
    
    
    if ~isempty(dir(ncfname))
        
        
        
        fprintf('Adding 2012 data  \n');
        
        [x_data, y_data, z_data ] = read_geomag_netcdf(ncfname,0,a);
        
        
        
        subplot(311);
        
        plot(time_array_min, x_data,'r');
        
        hold on;
        
        axis([time_array_min_2011(1) time_array_min(end) -inf inf]);
        
        datetick;
        
        ylabel('X');
        
        title(allstn(i,:));
        
        subplot(312);
        
        plot(time_array_min, y_data,'r');
        
        hold on;
        
        axis([time_array_min_2011(1) time_array_min(end) -inf inf]);
        
        
        
        datetick;
        
        ylabel('Y');
        
        subplot(313);
        
        plot(time_array_min, z_data,'r');
        
        hold on;
        
        axis([time_array_min_2011(1) time_array_min(end) -inf inf]);
        
        
        
        datetick;
        
        ylabel('Z');
        
        
        
        x_data(isnan(x_data)) = 99999.9;    %no data (this will get converted to 999999 by int32(99999.9*10)
        
        z_data(isnan(z_data)) = 99999.9;    %no data (this will get converted to 999999 by int32(99999.9*10)
        
        y_data(isnan(y_data)) = 99999.9;    %no data (this will get converted to 999999 by int32(99999.9*10)
        
        
        
        where_to_store = (datenum(2008,1,1) - datenum(1995,1,1)) * 24 * 60;
        
        length_of_record = a;
        
        
        
        netcdf.putVar(ncid_1, X_ID, where_to_store, length_of_record,  int32 (x_data .* 10) );
        
        netcdf.putVar(ncid_1, Y_ID, where_to_store, length_of_record,  int32 (y_data .* 10) );
        
        netcdf.putVar(ncid_1, Z_ID, where_to_store, length_of_record,  int32 (z_data .* 10) );
        
        
        
    end;
    
    
    
    netcdf.close(ncid_1);
    
    
    
    saveas(gcf,['/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/plots/' allstn(i,:)],'png');
    
    
    
    close all;
    
    
    
    
    
    
    
end;





%% create spline removed







S1 = dir('/nfs/satmag_work/mnair/projects/obs_mag_data/update_2011/absolute/*.nc');


S2 = dir('/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/work/combinednc/*.nc');


problem_stn = ['abg';'ebr';'lvv';'ppt';'thy';'vss'];



start_date = datenum(1995,1,1);

a = (datenum(2012,12,31) - datenum(1995,1,1) + 1) * 24 * 60;

time_array_min = (datenum(1995,1,1,0,0,0): (1/(24*60)): datenum(2012,12,31,23,59,0))';

basedirectory = '/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/absolute/';







% director for netCDF files

localdirectory = '/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/secular_removed/';



S1 = dir([basedirectory '*.nc']);





if exist('filestring','var')
    
    clear filestring;
    
end;






for i= 1:length(S1),
    
    filestring(i,:) = upper(S1(i).name(1:3));
    
end;


% Feb 5, 2012. For adding the quasi-dfinitive data from USGS, I use the
% following lines

filestring = upper(['brw';'bou';'cmo';'ded';'frd';'frn';'gua';'hon';'new';'sjg';'shu';'sit';'bsl';'tuc']);




[master_a,master_b] = xlsread('/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/work/Observatory_Master_Metadata.xlsx');



master_b(1,:) = [];



for i= 1:length(master_b),
    
    master_stn(i,:) = upper(cell2str(master_b(i,1)));
    
end;















for i = 1:length(filestring), % 55 GZH
    
    
    
    
    
    %create new netCDF file for writing the data
    
    
    
    [stn , index_a ,index_b ] = intersect(master_stn,filestring(i,:),'rows');
    
    
    
    metadatacellarray = master_b(index_a ,:);
    
    
    
    lat = master_a(index_a ,1);
    
    
    
    long = master_a(index_a ,2);
    
    
    
    alt = master_a(index_a ,3);
    
    
    
    fprintf('Creating %s ',filestring(i,:));
    
    
    
    [ncid_1, X_ID_1, Y_ID_1, Z_ID_1] = Create_Geomag_netCDF(localdirectory, lat, long, alt, metadatacellarray);
    
    
    
    
    
    %read absolute data
    
    
    
    ncfname = [basedirectory filestring(i,:) '_1995_2012.nc'];
    
    
    
    if ~isempty(dir(ncfname))
        
        
        
        fprintf('Reading Abs data  \n');
        
       
        % read the absolute data 
        [x_data, y_data, z_data , dummy, dummy, dummy, dummy] = read_geomag_netcdf(ncfname, 0, a, 0);
        
       
        
        
        
        
        
        L = isnan(x_data);
        
        data_array = x_data(~L);
        
        time_array = time_array_min(~L);
        
        
        
        b1 = min(time_array):365:max(time_array)+10;
        
        if length(b1) > 1,
            
            sp = spline(b1,data_array(1:60:end)'/spline(b1,eye(length(b1)),time_array(1:60:end)'));
            
            v = ppval(time_array_min,sp);
            
        else,%data length <=1 year
            
            sp = robustfit(time_array,data_array);
            
            v= sp(1) + sp(2) * time_array_min;
            
        end;
        
        
        
        % write sec removed x values
        
        
        
        subplot(311);
        
        plot(time_array_min,x_data);
        
        hold on;
        
        plot(time_array_min(~L),v(~L),'r');
        
        ylabel('X');
        
        title(filestring(i,:));
        
        
        
        axis([datenum(1995,1,1) datenum(2013,1,1) -inf inf]);
        
        datetick('x','keeplimits');
        
        x_data = x_data - v;
        
        
        
        x_data(isnan(x_data)) = 99999.9;    %no data (this will get converted to 999999 by int32(99999.9*10)
        
        
        
        
        
        netcdf.putVar(ncid_1, X_ID_1, 0, a,  int32 (x_data .* 10) );
        
        
        
        clear x_data v data_array time_array L;
        
        %y
        
        L = isnan(y_data);
        
        
        
        data_array = y_data(~L);
        
        time_array = time_array_min(~L);
        
        
        
        b1 = min(time_array):365:max(time_array)+10;
        
        if length(b1) > 1,
            
            sp=spline(b1,data_array(1:60:end)'/spline(b1,eye(length(b1)),time_array(1:60:end)'));
            
            v=ppval(time_array_min,sp);
            
        else,%data length <=1 year
            
            sp=robustfit(time_array,data_array);
            
            v= sp(1) + sp(2) * time_array_min;
            
        end;
        
        
        
        subplot(312);
        
        plot(time_array_min,y_data);
        
        hold on;
        
        plot(time_array_min(~L),v(~L),'r');
        
        ylabel('Y');
        
        
        
        
        
        axis([datenum(1995,1,1) datenum(2013,1,1) -inf inf]);
        
        datetick('x','keeplimits');
        
        % write sec removed y values
        
        
        
        y_data = y_data - v;
        
        
        
        y_data(isnan(y_data)) = 99999.9;    %no data (this will get converted to 999999 by int32(99999.9*10)
        
        
        
        netcdf.putVar(ncid_1, Y_ID_1, 0, a,  int32 (y_data .* 10) );
        
        
        
        clear y_data v data_array time_array L;
        
        
        
        %read abs z vlues
        
        
        
        
        
        L = isnan(z_data);
        
        data_array = z_data(~L);
        
        time_array = time_array_min(~L);
        
        
        
        %fit spline
        
        
        
        b1 = min(time_array):365:max(time_array)+10;
        
        if length(b1) > 1,
            
            sp=spline(b1,data_array(1:60:end)'/spline(b1,eye(length(b1)),time_array(1:60:end)'));
            
            v=ppval(time_array_min,sp);
            
        else,%data length <=1 year
            
            sp=robustfit(time_array,data_array);
            
            v= sp(1) + sp(2) * time_array_min;
            
        end;
        
        
        
        subplot(313);
        
        plot(time_array_min,z_data);
        
        hold on;
        
        plot(time_array_min(~L),v(~L),'r');
        
        ylabel('Z');
        
        
        
        
        
        axis([datenum(1995,1,1) datenum(2013,1,1) -inf inf]);
        
        datetick('x','keeplimits');
        
        % write sec removed z values
        
        
        
        z_data = z_data - v;
        
        
        
        z_data(isnan(z_data)) = 99999.9;    %no data (this will get converted to 999999 by int32(99999.9*10)
        
        
        
        
        
        netcdf.putVar(ncid_1, Z_ID_1, 0, a,  int32 (z_data .* 10) );
        
        
        
        
        
        
        
%         subplot(311);
%         
%         plot(time_array_min, x_data);
%         
%         axis([time_array_min(1) time_array_min(end) -inf inf]);
%         
%         
%         
%         datetick;
%         
%         ylabel('X');
%         
%         title(allstn(i,:));
%         
%         hold on;
%         
%         subplot(312);
%         
%         plot(time_array_min, y_data);
%         
%         axis([time_array_min(1) time_array_min(end) -inf inf]);
%         
%         
%         
%         datetick;
%         
%         ylabel('Y');
%         
%         hold on;
%         
%         subplot(313);
%         
%         plot(time_array_min, z_data);
%         
%         axis([time_array_min(1) time_array_min(end) -inf inf]);
%         
%         
%         
%         datetick;
%         
%         ylabel('Z');
%         
%         hold on;
%         
        
        
        
        
        
        
        
        
        
        
        
        
        
        
    end;
    
    
    
    netcdf.close(ncid_1);
    
    
    
    saveas(gcf,['/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/plots/' filestring(i,:) 'Sec_rem'],'png');
    
    
    
    close all;
    
    
    
    
    
    
    
end;



%% Add the USGS 2011 data to the database

usgsobs = ['brw';'bou';'cmo';'ded';'frd';'frn';'gua';'hon';'new';'sjg';'shu';'sit';'bsl';'tuc'];

% director for netCDF files

localdirectory = '/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/absolute/';


S1 = dir('/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/work/intermagnet/usgs_quasi_definitive/*.BIN');

datadirectory = '/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/work/intermagnet/usgs_quasi_definitive/';
start_date = datenum(1995,1,1);

for i= 1:length(S1),
    
    filestring(i,:) = upper(S1(i).name);
    
end;


[b, ia, ib] = unique(filestring(:,1:3),'rows');%Find the unique stations
a_2012 = (datenum(2012,12,31) - datenum(1995,1,1) + 1) * 24 * 60;



for i = 1:length(b),
    
    
    ncfname = [localdirectory b(i,:) '_1995_2012.nc'];
    
    [x_data, y_data, z_data , X_ID, Y_ID, Z_ID, ncid] = read_geomag_netcdf(ncfname, 0, a_2012, 1);
    
    
    
    L = ib == i; %select a station
    
    filesubset = filestring(L,:);
    
    for j = 1:sum(L),
        
        intermagfname = [datadirectory filesubset(j,:)];
        
        data_struct = INTERMAG(intermagfname, 'XYZF');
        
        
        
        
        
        if ~isempty(data_struct) && (max(abs(diff(data_struct.MINUTELY(:,1)))) < 1.1/(24*60)),
            
            
            
            data_struct.MINUTELY(isnan(data_struct.MINUTELY)) = 99999.9;    %no data (this will get converted to 999999 by int32(99999.9*10)
            
            fprintf('file name: %s\n',intermagfname);
            
            where_to_store = ( data_struct.MINUTELY(1,1) - start_date ) * 1440 ;
            
            length_of_record = length(data_struct.MINUTELY(:,1));
            
            
            
            %netcdf.putVar(ncid,varid,start,count,data)
            
            %multiplied by 10 to preserve the decimal points
            
            netcdf.putVar(ncid, X_ID, where_to_store, length_of_record,  int32 (data_struct.MINUTELY(:,2) .* 10) );
            
            netcdf.putVar(ncid, Y_ID, where_to_store, length_of_record,  int32 (data_struct.MINUTELY(:,3) .* 10) );
            
            netcdf.putVar(ncid, Z_ID, where_to_store, length_of_record,  int32 (data_struct.MINUTELY(:,4) .* 10) );
            
            fprintf('%s             OK \n',intermagfname);
            
        else
            
            fprintf('Error %s\n',intermagfname);
        end;
    end;
    
    netcdf.close(ncid);
    
end;


%% Read absolute data, plot data & spline, save the figure
% also write the data attributes to a file


S1 = dir('/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/absolute/*.nc');


start_date = datenum(1995,1,1);
a = (datenum(2012,12,31) - datenum(1995,1,1) + 1) * 24 * 60;
time_array_min = (datenum(1995,1,1,0,0,0): (1/(24*60)): datenum(2012,12,31,23,59,0))';
basedirectory = '/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/absolute/';
S1 = dir([basedirectory '*.nc']);

fid = fopen('/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/station_list_new.txt','wt');

if exist('filestring','var')
    clear filestring;
end;

x_data_av_matrix = nan([155,263]);
y_data_av_matrix = nan([155,263]);
z_data_av_matrix = nan([155,263]);

for i= 1:length(S1),
    filestring(i,:) = upper(S1(i).name(1:3));
end;

[master_a,master_b] = xlsread('/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/work/Observatory_Master_Metadata.xlsx');
master_b(1,:) = [];
for i= 1:length(master_b),
    
    master_stn(i,:) = upper(cell2mat(master_b(i,1)));
    
end;

for i = 1:length(filestring), % 55 GZH
    
    
    %create new netCDF file for writing the data
    
    [stn , index_a ,index_b ] = intersect(master_stn,filestring(i,:),'rows');
    
    metadatacellarray = master_b(index_a ,:);
    
    lat = master_a(index_a ,1);
    
    long = master_a(index_a ,2);
    
    alt = master_a(index_a ,3);
    
 
    
    
    %read absolute data
    
    ncfname = [basedirectory filestring(i,:) '_1995_2012.nc'];
    
    if ~isempty(dir(ncfname))
        
        fprintf('Reading %s data  \n', filestring(i,:));
        % read the absolute data 
        [x_data, y_data, z_data , dummy, dummy, dummy, dummy] = read_geomag_netcdf(ncfname, 0, a, 0);
      
        
        L = isnan(x_data);
        
%         data_array = x_data(~L);
%         
%         time_array = time_array_min(~L);
%         
%         b1 = min(time_array):365:max(time_array)+10;
%         
%         if length(b1) > 1,
%             
%             sp = spline(b1,data_array(1:60:end)'/spline(b1,eye(length(b1)),time_array(1:60:end)'));
%             
%             v = ppval(time_array_min,sp);
%             
%         else,%data length <=1 year
%             
%             sp = robustfit(time_array,data_array);
%             
%             v= sp(1) + sp(2) * time_array_min;
%             
%         end;
        
        max_x_data = sum(~L);
        x_data_av_matrix(i,:) = any(reshape(~L,[1440*25,263]),1);
        
        % write sec removed x values
        
        
        
%         subplot(311);
%         
%         plot(time_array_min,x_data);
%         
%         hold on;
%         
%         plot(time_array_min(~L),v(~L),'r');
%         
%         ylabel('X');
%         
%         title(filestring(i,:));
%         
%         
%         
%         axis([datenum(1995,1,1) datenum(2013,1,1) -inf inf]);
%         
%         datetick('x','keeplimits');
        
%         x_data = x_data - v;
        
        
        
%         x_data(isnan(x_data)) = 99999.9;    %no data (this will get converted to 999999 by int32(99999.9*10)
%         
%         clear x_data v data_array time_array L;
        
        %y
        
        L = isnan(y_data);
        
        
        
%         data_array = y_data(~L);
%         
%         time_array = time_array_min(~L);
%         
%         
%         
%         b1 = min(time_array):365:max(time_array)+10;
%         
%         if length(b1) > 1,
%             
%             sp=spline(b1,data_array(1:60:end)'/spline(b1,eye(length(b1)),time_array(1:60:end)'));
%             
%             v=ppval(time_array_min,sp);
%             
%         else,%data length <=1 year
%             
%             sp=robustfit(time_array,data_array);
%             
%             v= sp(1) + sp(2) * time_array_min;
%             
%         end;
%         
        max_y_data = sum(~L);
        y_data_av_matrix(i,:) = any(reshape(~L,[1440*25,263]),1);
   
        
%         subplot(312);
%         
%         plot(time_array_min,y_data);
%         
%         hold on;
%         
%         plot(time_array_min(~L),v(~L),'r');
%         
%         ylabel('Y');
%         
%         
%         
%         
%         
%         axis([datenum(1995,1,1) datenum(2013,1,1) -inf inf]);
%         
%         datetick('x','keeplimits');
%         
%         % write sec removed y values
%         
%         
%         
%         y_data = y_data - v;
%         
%         
%         
%         y_data(isnan(y_data)) = 99999.9;    %no data (this will get converted to 999999 by int32(99999.9*10)
%         
%         
%         
%     
%         
%         
%         clear y_data v data_array time_array L;
        
        
        
        %read abs z vlues
        
        
        
        
        
        L = isnan(z_data);
        
%         data_array = z_data(~L);
%         
%         time_array = time_array_min(~L);
        
        
        
        %fit spline
        
        
        
%         b1 = min(time_array):365:max(time_array)+10;
%         
%         if length(b1) > 1,
%             
%             sp=spline(b1,data_array(1:60:end)'/spline(b1,eye(length(b1)),time_array(1:60:end)'));
%             
%             v=ppval(time_array_min,sp);
%             
%         else,%data length <=1 year
%             
%             sp=robustfit(time_array,data_array);
%             
%             v= sp(1) + sp(2) * time_array_min;
%             
%         end;
%         
             max_z_data = sum(~L);
           z_data_av_matrix(i,:) = any(reshape(~L,[1440*25,263]),1);

        
%         subplot(313);
%         
%         plot(time_array_min,z_data);
%         
%         hold on;
%         
%         plot(time_array_min(~L),v(~L),'r');
%         
%         ylabel('Z');
%         
%         
%         
%         
%         
%         axis([datenum(1995,1,1) datenum(2013,1,1) -inf inf]);
%         
%         datetick('x','keeplimits');
        
        % write sec removed z values
        
        
        
%         z_data = z_data - v;
%         
%         
%         
%         z_data(isnan(z_data)) = 99999.9;    %no data (this will get converted to 999999 by int32(99999.9*10)
        
        
    end;
    
    
%     saveas(gcf,['/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/plots/' filestring(i,:) 'Sec_rem'],'pdf');
    
    fprintf(fid,'%4s %9.4f % 9.4f %6.1f %d\n',filestring(i,:), lat, long, alt, min([max_x_data, max_y_data, max_z_data]));
    
    close all;
    
    
    
    
    
    
    
end;

fclose (fid);

%% The following script was used to make a data availability plot for the report

 load /nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/data_availability x_data_av_matrix y_data_av_matrix z_data_av_matrix filestring
 
 imagesc(x_data_av_matrix & y_data_av_matrix & z_data_av_matrix);
 
 % now plot the station names
 set(gca,'YTick',[]);
for i = 1:2:155,
text(-25,i+0.5, filestring(i,:),'FontSize',3);
end;

  for i = 2:2:155,
text(270,i+0.5, filestring(i,:),'FontSize',3);
end;

set(gca,'XTick',1:263/18:264);
set(gca,'XTickLabel',['1995';'1996';'1997';'1998';'1999';'2000';'2001';'2002';'2003';'2004';'2005';'2006';'2007';'2008';'2009';'2010';'2011';'2012';'2013']);
